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METHOD FOR IMAGING MULTIPHASE FLOW USING ELECTRICAL X> 
CAPACITANCE TOMOGRAPHY 

DESCRIPTION 
PRIORITY 

[0001] This application claims the benefit under 35 U.S.C. § 1 19(a)-(d) or § 365(a)- 

(c), of an application filed under the Patent Coorporation Treaty ("PCT"), no. 
PCT/MX2003/000067 filed August 22 2003, the entire contents of which are hereby 
incorporated by reference. 



TECHNICAL FIELD OF INVENTION 

[0002] This invention relates to a new image reconstruction method for imaging 

multiphase flows using electrical capacitance tomography (ECT), it will be used in 
tomographic multiphase flow meters to quantify the various fluids (such as gas, oil and 
water) produced by an oil well. It can also be used to image and monitor other multiphase 
flows and processes occurring in industry, such as: optimization of design and operation of 
separator tanks, optimization of design of fluidized catalytic cracking units (FCC) and bed 
fluidized systems, and optimization of design of bed reactors. 

BACKGROUND OF INVENTION 

[0003] In general terms, tomography is used for obtaining an image of a cross section 

of an object in a given plane. X-ray tomography was the first to be developed (in 1970s) and 
its use is now routine not only in medicine but in some industrial applications as well 
(internal inspection of mechanical components and flaw detection in materials, for example). 
[0004] Subsequently, a number of new tomographic methods aimed at industrial 

processes have emerged, collectively known as process tomography (Williams and Beck, 
1995). The aim of these methods, which started to develop in the mid 1980s, is to produce an 
image of the phase or component distribution in an industrial process using only external 
sensors and without causing any perturbation to it, as depicted in figure 1, which shows a 
system of process tomography formed by a tank or pipe (A), a data acquisition system (B), 
and a computer for image reconstruction (C). In other words, process tomography provides a 
way of 'looking' from the outside to inside the process, with no need for physical intrusion or 
alteration, and this represents a radically different approach to gathering structural 
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information on the process to a global scale, unlike the traditional methods based on local 
sampling on a certain number of points. The process may occur, for instance, in mixing or 
stirring vessels, fluidized bed reactors, separator tanks, or a pipeline carrying multiphase 
flow. 

[0005] There is a whole range of principles and techniques that can be exploited in 

process tomography, including electrical methods based on impedance measurement, 
ultrasound, magnetic resonance, optical methods and those based on ionizing radiation (X- 
and gamma-rays). Generally speaking, ionizing radiation methods produce images with the 
highest definition, but are relatively slow. On the other hand, electrical methods yield low- 
resolution images but are much faster, robust and relatively inexpensive. 
[0006] In particular with regard to electrical impedance tomography, or electrical 

tomography for short, there has been a very noticeable progress in the last few years. This 
type of tomography has two main modalities: capacitance tomography and resistance 
tomography. In a capacitance tomography system (Beck et aL 9 1997; Gamio, 1997; 
Plaskowski et aL, 1995), (as depicted in figure 2, formed by a sensor (A), a data acquisition 
system (B) and an image reconstruction computer (C)), normally used with mixtures where 
the continuous phase is non-conducting, the sensor employed, as depicted in figures 3 and 4, 
is made of a circular array of electrodes (B) distributed around the cross-section to be 
examined, and the capacitance between all the different electrode-pair combinations is 
measured. With the help of a computer and a suitable image reconstruction algorithm, this 
information is used to create a map (an image) showing the variation of the dielectric 
constant (or relative permittivity) inside the sensor area, thus providing an indication of the 
physical distribution of the various components of the mixture. Relative to figures 3 and 4, 
the electrodes (B) can be located on the outside of a non-conducting pipe (A), in order to 
simplify sensor construction and avoid direct contact with the process fluids (E). A second 
external grounded metallic pipe (C) serves as an electric screen and to provide mechanical 
resistance. The sensor also has two cylindric guard electrodes (D) that are grounded as well. 
[0007] The British patent GB 2 214 640, issued September 6, 1989, describes an 

electrical capacitance tomography system that employs a linear back-projection (LBP) 
algorithm as the image reconstruction method. However, such reconstruction method 
produces images of a relatively low definition. 

[0008] Resistance tomography, on the other hand, is aimed at mixtures where the 

continuous phase is a conductor of electricity (Plaskowski et aL, 1995; Williams y Beck, 
1995). In this case, the electrodes are installed flush with the inside surface of the pipe (or 
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vessel) wall and in direct contact with the fluids. A number of different excitation current 
patterns are applied and the resultant voltages are measured. They are then used to construct a 
map of the conductivity distribution inside the sensor, which reflects the physical distribution 
of the mixture components. 

[0009] In principle, electrical capacitance tomography (ECT) has important 

applications in multiphase flow measurement, particularly gas-oil two-phase flow, which 
often occurs in many oil wells. The traditional way to quantify the various fluids produced by 
an oil well is to separate the mixture by gravity in large tanks, prior to measuring each 
component separately using conventional single-phase flow meters. In the last decades, 
multiphase flow meters have appeared which allow the quantification of the produced fluids 
without the need to separate the mixture (Thorn et al. 9 1997). However, the multiphase meters 
currently available suffer from an unwanted sensibility to changes in the flow regime, unless 
they are equipped with flow mixing or conditioning devices that introduce permanent 
pressure losses (which ultimately translate into energy loss). This limitation should be 
avoided through ECT as the flow regime could be determined and used to compensate the 
response of conventional multiphase meters, or, alternatively, it is possible to design a new 
type of tomographic multiphase meter, based on analyzing series of ECT images from two 
slightly separated cross-sections of the pipe (Hammer et al. 9 1997; Plaskowski et aL, 1995). 
Additionally, ECT has potential applications to imaging, monitoring and controlling 
numerous industrial multiphase processes. 

[0010] However, so far the main limiting factor to the practical application of ECT 

has been the lack of fidelity or accuracy of the images obtained using the available image 
reconstruction methods, giving rise to a need of improved methods as the one introduced in 
this invention (Yang y Peng, 2003). Simple direct methods like linear back-projection (LBP) 
yield relatively poor images that only provide a qualitative indication of the component 
distribution inside the sensor. One said method is described in British patent GB 2214640, 
issued September 6, 1989. 

[0011] On the other hand, more sophisticated methods, based on iterative local 

optimization techniques, generally require one or more regularization parameters whose 
optimal value depends precisely on the (unknown) image to be reconstructed, apart from the 
fact that the regularization employed has the effect of smoothing the image contours, making 
it more diffuse. One said method is described in British patent GB 2329476 A (issued March 
24, 1999). 
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[0012] Thus, there is an urgent need of better and more accurate image reconstruction 

methods. As an example of this, a patent was issued recently covering an image 
reconstruction method based on the application of artificial neural networks (American patent 
US 6,577,700 Bl, issued July 10, 2003). 

[0013] The present invention describes new image reconstruction methods based on 

simulated annealing and genetic algorithms. According to the laws of physics, electrostatics 
in particular, the ECT sensor (as can be observed in figure 4) can be considered as a special 
case of a system of charged conductors separated by a dielectric medium, the theory of which 
was first developed by J. C. Maxwell (1873). In our particular case, the sensor electrodes (B) 
act as the charged conductors while the insulating pipe end (A) the sensor contents acts as the 
dielectric medium. For an n electrode sensor, the induced electrode charges q ( and the 
electrode potentials v, are related by the following set of linear equations 

q x = c u v x + c 12 v 2 + + 
q 2 = c 2I v, + c 22 v 2 + + 

<l n = c ml v t + c n2 v 2 + + 

[0014] where c« are the self-capacitance coefficients (or just self-capacitances for 

short) of electrode z, while the others, , with i ^ j 9 are the mutual capacitance coefficients 
(or just mutual capacitances) between electrodes i and j. Put in matrix form, equation (1) 
becomes 

3* 

or 

q = Cv (3) 

[0015] Since the capacitances have the property of reciprocity, i.e., eg = cji , there are 

only m = \n{n-\) independent mutual capacitances, corresponding to lower (or upper) 
triangular matrix of C, and corresponding also to each one of the m different electrode- pairs 
that can be formed in the sensor. 
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[0016] The value of the mutual capacitances is a complex non-linear function of the 

conductor system geometry, and of the spatial distribution of the dielectric constant or 
relative permittivity (hereafter called just 'permittivity 1 ) of the dielectric medium. In the case 
of the ECT sensor, the geometry of the electrodes, that of the pipe, and the value of the 
dielectric constant of the latter, are all fixed. Therefore, it can be said that the mutual 
capacitances are a function only of the spatial distribution of the dielectric constant inside the 
sensor, s(x 9 y). The problem of calculating the mutual capacitances corresponding to a specific 
permittivity distribution inside the sensor is referred to as the forward problem. 
[0017] The use of the cylindrical end guards (figure 3) and the assumption that the 

phase (and thus the permittivity) distribution does not change too much in the axial direction, 
allows the sensor to be represented by a two-dimensional (2-D) model (Xie et al. 9 1989). 
Unless otherwise stated explicitly, in what follows a 2-D model of the sensor will be used. 
Therefore, the electric charges q ( and the capacitances c, 7 y c,y that appear in equations (1) and 
(2) should be considered quantities per unit length of the electrodes in the axial direction. A 
tilde (i.e., q i , c H and c 0 ) shall be used to denote the total quantities that result from 

considering the actual length of the electrodes. The previous variables are related between 
them through the electrode length L 9 according to 

<7, = J . c. =2l and c u = Y ( 4 ) 

[0018] If the interior of the 2-D sensor is divided into p equal-area regions (or 

'pixels') where the permittivity is considered constant, then the discrete version of the forward 
problem is 





c \ 




'/.(e)' 






c = 




= f(B) = 


./-(e). 


with e = 





[0019] where c is the vector of mutual capacitances (per unit length),/ are non-linear 

functions not known explicitly and s is the vector of permittivities corresponding to the p 
regions or pixels within the sensing zone. 

[0020] Applying Gauss's Law, the mutual capacitances per unit axial electrode 

length can be calculated as 
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[0021] where s Q is a physical constant called the permittivity of free space, equal to 

8.854 x 10" 12 farads per meter, Si is a closed curve surrounding electrode i 9 dl is a normal 
vector representing an element of the curve T, > dl is an element of length of that curve, the 
symbol 1 • ' represents the scalar product of two vectors, and ft is the electrostatic potential 
produced in the sensor when applying a voltage of V volts to electrode j (which is called 
source or excitation electrode) and 0 volts to all others (called detection electrodes). 
[0022] The potential ft is determined by the solution of the following partial 

differential equation 

V€{x,y)V(f> J =0 (7) 

[0023] subject to the boundary conditions (a) <fjf = V volts on the source electrode and 

(b) ft = 0 on the detection electrodes and the outer screen. In general, equation (7) does not 
have an analytic solution and must be solved numerically. 

[0024] The problem of estimating what is the spatial permittivity distribution 

inside the sensor that corresponds to a specific set of mutual capacitance values is referred to 
as the inverse problem, and is the problem that image reconstruction methods must address 
and solve. Normally, the permittivity estimation is made in a discrete way, representing it as a 
vector 8 like the one in equation (5), which must be calculated from a vector of observed 
mutual capacitances c, obtained using of a suitable measurement apparatus. 
[0025] In order to solve the inverse problem, most ECT systems employ the 

linear back-projection algorithm (LBP) (Plaskowski et al., 1995; Yang y Peng, 2003; Xie et 
aL, 1989, 1992), which is described next. As a first step, a sensitivity map must be calculated 
for each one of the m = kn(n-\) possible electrode pairs, given by 

J,0fc) = C ' W ~ C '<^> for / = 1, . . , n (8) 

C /(>//) C i(emp) 

where k is the pixel number (from 1 to /?), c ( {k) is the capacitance measured with electrode 
pair i when the area of pixel k is full of a high-permittivity material while the rest of the 
sensor is full of a low-permittivity material, whereas ctyw/) and c^ emp ) are the capacitances for 
electrode pair i when the sensor is full of high- and low-permittivity material, respectively. 
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Generally, these sensitivity maps are calculated by solving numerically equation (7) and 
applying equation (6). 

[0026] Having determined the sensitivity maps, they can be used to obtain a 

permittivity image from any vector ofw = %n{n-\) measured mutual capacitances, c, 
corresponding to some unknown material distributions inside the sensor. For this, the 
measured capacitance readings must first be normalized according to 

= C ' ~ C/ (ewp) (9) 

C i(full) ~ C i(emp) 

[0027] where A, is the normalized capacitance for electrode pair i and c, is the actual 

capacitance measured with that electrode pair. 

[0028] The basic LBP formula calculates a 'grey level 1 g{k) for each pixel as 

Z>u(*) 

*(*) = for *= i, . ■ ,p (io) 

1=1 

In principle, this grey level is supposed to be linearly related to the permittivity, with g=l 
and g = 0 corresponding to the permittivities of the high- and low-permittivity materials, 
respectively. LBP is based on making a linear approximation to a problem that, as already 
mentioned, is essentially non-linear (Gamio and Ortiz- Aleman, 2003). Therefore, this image 
reconstruction method causes considerable errors, which are particularly grave if there are 
large differences in permittivity in the image. 

[0029] So far, the main alternative to LBP has been the use of iterative methods that 

seek to minimize some objective function, employing local optimization techniques like the 
regularized Newton-Raphson method or other similar approaches (Yang and Peng, 2003). As 
an example of these methods, there is the one used in the EIT2D software package 
(Vauhkonen et al. 9 2001), developed by researchers from Finland and the United Kingdom. 
Their method is based on minimizing, with respect to e, the following functional 

\\c meas -c calc \\ 2 + ^IILeU 2 (11) 

where a is a regularization parameter, L is a 'regularization' matrix containing some type of 
ci-prion smoothness information about c, and c ca i c — f(e) is the vector of n calculated mutual 
capacitance values for a given permittivity distribution inside the sensor. Starting with an 
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initial guess e 0 , the minimization is carried out by the following iterative procedure (basically 
a Newton-type method with Tikhonov regularization) 

«»♦, = e* + [J/J* + a 2 L T L]-' {J/ [c meas - f(e,) ] - a 2 L T Le, } (12) 
where J* is the so-called Jacobian matrix of the partial derivatives of f(e), evaluated at 



J, = J(e,) = 



K J 



( f=l...m, 7=l.../7) (13) 



[0030] However, these image reconstruction methods have the problem that they 

require, for their correct operation, one or more regularization parameters whose right value 
is strongly dependent, precisely, on the image that one wishes to reconstruct, implying that 
one would need to know beforehand the solution to the problem. Moreover, these methods 
produce distorted images, because the regularization has an excessive smoothing effect on the 
obtained permittivity. If the regularization is too strong the smoothing effect will occur, and if 
it is too weak the method can become unstable and/or not converge to the desired solution. 
[0031] These local optimization algorithms, during their search, explore only a 

relatively small sector of the solution domain, restricted to the vicinity of the initial guess. If 
the optimal solution of the problem, i.e., the absolute minimum of the objective function, is 
located far from the initial guess, it will hardly be reached due to the presence of relative 
minima in their way, places where these methods can become trapped. The most used 
methods in this category are least-squares linear inversion and techniques that utilize the 
gradient of the objective function, like the steepest-descent and the conjugate-gradient 
methods. In general, local search methods exploit the (scarce) information derived from the 
comparison of a small number of models (solutions), thus avoiding an extensive search in the 
whole model space (Sambridge y Drijkoningen, 1992). 

[0032] Global optimization methods explore the whole solution domain during the 

inversion process. They carry out an extensive scan within the model space. In this way, 
despite the existence of partial solutions to the problem, there is a greater possibility that the 
final solution corresponds to the best fit between the observed data and the synthetic ones. 
This type of methods, contrary to local techniques, does not require the information provided 
by the derivatives of the objective function, because in this case the problem does not need to 
be linearized. Global optimization algorithms use stochastic criteria in order to 
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simultaneously explore all the solution space in search of the optimal model. The best known 
of the global methods is Monte Carlo, which performs a purely random and unbiased search. 
In other words, when generating each new model, it does not take advantage of the 
information obtained from the previously evaluated models (Gallagher et aL, 1991). The 
unguided randomness is the most characteristic feature of this method, which distinguishes it 
from the rest of the global methods. Among the global optimization techniques, there are also 
the method of simulated annealing and genetic algorithms. Both were conceived as analogies 
of optimization systems occurring in nature. Genetic algorithms emulate the mechanisms of 
biological evolution while simulated annealing is based on thermodynamics. Both methods 
are inherently non-linear and, therefore, lend themselves naturally to their application in 
capacitance tomography, a non-linear problem. 

SPECIFICATION 

[0033] The present invention refers to an image reconstruction method for imaging 

multiphase flows using electrical capacitance tomography (ECT), based on the use of 
heuristic non-linear techniques of global optimization, particularly the simulated annealing 
method and genetic algorithms. It considers a circular array of rectangular contiguous 
metallic electrodes placed around the outer wall of a pipe made of an electrically insulating 
material, forming a sensor. Trough the sensor flows a mixture of fluids in the form of a 
multiphase flow, whose spatial distribution inside the sensor is to be determined. Data are 
obtained by performing mutual capacitance measurements between all possible electrode 
pairs. Said data depend on the fluid distribution inside the pipe. 

[0034] A matter of this invention is to provide a method to process said data in order 

to reconstruct an image of the spatial distribution of the phases or components of the 
multiphase mixture that flows through the sensor, using the method of simulated annealing 
and/or genetic algorithms. 

[0035] Another matter of this invention is that it can be used in tomographic 

multiphase flow meters to quantify the various fluids (such as gas, oil and water) produced by 
an oil well. 

[0036] Another matter of the present invention is that it can also be used to image and 

monitor other multiphase flows and processes occurring in industry, such as optimization of 
design and operation of separator tanks, optimization of design of fluidized catalytic cracking 
units (FCC) and bed fluidized systems, and optimization of design of bed reactors. 
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BRIEF DESCRIPTION OF THE DRAWINGS 

[0065] Figure 1 represents schematically a process tomography system employed in 

the present invention. 

[0066] Figure 2 represents schematically an electrical capacitance tomography system 

employed in the present invention. 

[0067] 

[0068] Figure 3 shows a drawing of the sensor employed in capacitance tomography 

used in the present invention. 

[0069] Figure 4 shows a cross-section of the sensor at the zone of the measuring 

electrodes employed in the present invention. 

[0070] Figure 5 is a schematic drawing illustrating the method of simulated annealing 

employed in the present invention. 

[0071] Figure 6 shows the flow diagram of one possible implementation of the 

method of simulated annealing employed in the present invention. 

[0072] Figure 7 is a schematic drawing illustrating the method of genetic algorithms 

employed in the present invention. 

[0073] Figure 8 shows schematically the process of crossover used in genetic 

algorithms employed in the present invention. 

[0074] Figure 9 shows schematically the process of mutations used in genetic 

algorithms employed in the present invention. 
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[0075] Figure 10 shows the flow diagram of one possible implementation of the 

method of genetic algorithms employed in the present invention. 

[0076] Figure 1 1 shows the results obtained by reconstructing images of a simulated 

gas-oil annular flow, using the simulated annealing method in the present invention. 
[0077] Figure 12 shows the results obtained by reconstructing images of a simulated 

gas-oil stratified flow, using the simulated annealing method in the present invention. 
[0078] Figure 13 shows the results obtained by reconstructing images of a simulated 

gas-oil bubbly flow, using the simulated annealing method in the present invention. 

DETAILED DESCRIPTION 

[0079] The new image reconstruction procedures of the present invention are based 

on heuristic global optimization methods, specifically simulated annealing and genetic 
algorithms. 

[0080] A plurality of n metallic electrodes is placed around the periphery of a region 

to be imaged and electrical measurements of capacitance or resistance are collected between 
them. That is to say, m = \n(n-\) measurements (or data) are obtained. It is preferred that the 
data is capacitance data but it can also be resistance data. As shown in figure 3, the electrodes 
(B) can be rectangular in shape and be placed on the outer wall of a pipe (A) made of 
electrically non-conducting material, thus forming a sensor, which contains a multiphase or 
multicomponent flow (E). The aim is to use said measurement data to reconstruct an image of 
the spatial distribution of the dielectric constant (or permittivity) within the imaged region. 
Said permittivity distribution reflects the distribution of the substances filling the interior of 
the sensor, where the imaging area is. This area shall be considered divided into p equal parts. 
[0081] Method of Simulated Annealing 

[0082] The simulated annealing method is based on an analogy with the 

thermodynamic process of crystallization. A mineral fluid that cools slowly until it reaches a 
low energy state, gives rise to the formation of well defined crystals. If, on the contrary, the 
substance leaves its thermal equilibrium state with a sudden or partial cooling, the resulting 
crystal will have many defects, or the substance may even form a 'glass 1 , characterized by its 
meta-stable molecular disorder. This concept is used in the context of optimization methods 
to recognize potentially useful models or configurations. 

[0083] The atoms of each molecular configuration are equivalent to the model 

parameter in the inverse problem (i.e., the permittivity of the various image pixels). The 
system energy for such configuration is related to the cost (or misfit) function associated with 
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the set of parameters involved in the model. In our case, the system energy is associated with 
the following L2 norm 

m 

[0084] E = L 2 = (1 = 1, . . , m) (14) 

where c(i) meas are the m measured capacitances and c(i) ca i c are the ones calculated by solving 
the forward problem for a given permittivity distribution e. From an initial permittivity 
distribution, the method generates a range of configurations or parameter combinations 
considering a certain temperature Tfor the process. For this purpose the Metropolis et al. 
(1953) criterion is employed, which consists in changing a parameter, in each iteration, by a 
small random amount. This shift causes a change AE in the system's total energy. If AE is less 
than or equal to zero, the change in the parameter is accepted and the resulting configuration 
is considered as the new current configuration. When there is an increase in the system 
energy (AE is greater than zero), the probability of acceptance or rejection for the parameter 
change is determined as 

[0085] P(AE) = e-* E/T (15) 

[0086] In order to decide whether or not a change that produces an increase in the 

system energy is accepted, a random number between cero and one is chosen, which is then 
compared with the value of the probability corresponding to AE. If said random number is 
smaller, the parameter shift is accepted and the new configuration is considered as the current 
(updated) one. If said random number is greater, the parameter shift is not accepted and the 
configuration that existed before the shift is maintained. Repeating this procedure 
continuously, the thermal movement of the atoms of a system in thermal equilibrium (at a 
fixed temperature T) is simulated. I order to reach the system's base state, that is to say, the 
state of lowest energy and highest order, the temperature must be reduced very slowly, 
simulating a quasi -static process. This means that, during the cooling, the system must 
experience a series of states infinitesimally separated from the state of thermal equilibrium. 
[0087] The method of simulated annealing has three basic components (Vasudevan, 

1991): an energy (or cost, or misfit) function, an order function (the Metropolis criterion), 
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and a parameter that controls the system temperature. The process consists of three nested 
cycles. Figure 5 shows a diagram that illustrates the method of the present invention. 

The external cycle (3) regulates the system temperature. Every time a cycle is completed, the 
temperature decreases as it is multiplied by a factor R T thai is normally very close to one 
(0 < R T < 1). In this way the desired slow and gradual cooling is carried out. The intermediate 
cycle (2) updates the values, independent of each other, of a series of constants K ( associated 
with each parameter. Said constants determine the maximum change that each parameter may 
experience when it is perturbed in the innermost cycle (1). The value of said constants 
depends on the number of times that the current model has been accepted (according to the 
Metropolis criterion) at the end of every sequence of internal cycles (1). In the internal cycle 

(1) the parameter values are perturbed using the factors Ki , defined in the intermediate cycle 

(2) . The perturbation is done multiplying each parameter by the product of its corresponding 
Ki times a randomly chosen number between minus one and one. After this, the synthetic 
response of the current model is calculated and the change in the system's energy associated 
with the new parameter configuration is evaluated. Said energy change corresponds to the 
misfit between the synthetic data curve and the observed or measured one. If the misfit 
decreases, then the new configuration will be accepted as the current one and in turn 
perturbed in the same way. If, on the contrary, if the random perturbation causes an increase 
in the misfit, associated with an increment in the energy E, then that configuration is assigned 
a probability of acceptance according to the Metropolis criterion. 

[0088] The cycles (1), (2) and (3) of figure 5 are repeated, while the temperature of 

the process decreases progressively. As the temperature diminishes, the parameter variations 
are smaller and smaller. In this way, the search in the solutions domain tends to confine itself 
towards the models associated with the absolute minimum of the misfit function E. The end 
result is a set of values for the parameters (i.e., the permittivity in the various pixels that make 
an image) whose synthetic response reproduces the observed (capacitance) data, with a 
sufficiently small error. 

[0089] As an example, one possible specific embodiment of, albeit not the only one, 

of the method of simulated annealing, is presented in figure 6. The method starts in block (1), 
with a series of initial values for the temperature, the perturbation constants, the permittivities 
and the cost or misfit function (T Q , AT/( 0) y s^ 0 ) (con i = 1, ...,/?) y E Q ), and initializing the 
external cycle counter, denoted by k. Next, in block (2) the counters for the internal and 
intermediate cycles, i and j 9 are initialized, and the internal cycle starts. In it, one by one the 
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parameters (or permittivities) are subject to a random perturbation in block (3). Still in block 
(3), each time a parameter is perturbed, the forward problem is solved and the cost or misfit 
function, E, is calculated applying equation (14). If there was a decrease of E with respect to 
the previous evaluation of it, the perturbed parameter value is accepted as the new current 
value, the internal cycle parameter counter i is increased by one, and the method proceeds to 
perturb the next parameter (if there is one). If, on the contrary, there was an increase in E, 
then the Metropolis criterion is applied, in block (4), in order to decide whether or not the 
perturbed parameter value is accepted as the new current one. If, according to said criterion, 
the new value is accepted, then the internal cycle counter i is incremented by one and the 
method proceeds to perturb the next parameter (if there is one). If, according to the 
metropolis criterion, the perturbed value is not accepted, then the counter i is not incremented 
and the method proceeds to perturb once again the same parameter. Once all parameters Si 
have been updated, in block (5) the value of the constants AT, (which determine how the 
parameters are perturbed within the internal cycle) is adjusted, and the intermediate cycle 
counter j is increased by one. Nk determines how many times the internal cycle will be 
repeated without reducing the temperature. In other words, the intermediate cycle consists in 
the repetition of the internal cycle N K times, but with different values for AT, . This is done so 
in order to prevent the temperature from descending too fast, which can have negative effects 
in some simulated annealing applications. However, in the present invention that does not 
happen, and the value of N K is taken as one. At the end of the intermediate cycle, in block (6) 
the temperature is reduced as indicated in previous paragraphs and the external cycle counter 
k is increased by one. The whole previous procedure is repeated until k reaches the iteration 
limit Nt, or before that if a sufficiently low value for the cost function E is obtained. 
[0090] The present invention also refers to a image reconstruction method based on 

genetic algorithms, that is described in the next section. 
[0091] Method of Genetic Algorithms 

[0092] Genetic algorithms, originally proposed by Holland (1975), represent an 

evolution of the Monte Carlo method for strongly non-linear problems. The search for the 
optimum model is carried out exploring simultaneously the whole solution space, employing 
a probabilistic transition rule to guide the search. The process starts from a set of randomly 
chosen models. 

[0093] The parameters of each model are transformed into binary code in 

order to form chains called chromosomes, which are then subject to natural and genetic 
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selection criteria. The processes of selection, crossover and mutation update the population of 
models, originating a new generation of chromosomes, emulating the way in which biological 
evolve to produce organisms better adapted to the environment. The whole process is 
repeated until the measure of the misfit function approaches the maximum fit for all the 
population. 

[0094] The flow diagram of figure 7 summarizes the process utilized to apply 

an inversion scheme based on genetic algorithms, similar to the one described by Rodriguez- 
Zuniga et al. (1996) and Ortiz-Aleman et aL (2002, 2003) for the inversion of underground 
elastic parameters from surface inclination data and for the inversion of aeromagnetic data, 
respectively. The basic steps of the process are briefly described next, with reference to figure 
7. 



[0095] Discretization 

[0096] In the block (1) of figure 7, the parameters are represented by the vector of 

unknowns 6, which represents the unknown permittivity distribution. The cost function, 
which determines the fit between the observed data and the synthetic response of a given 
model, is denoted by E(e). The parameter coding is done taking into account the required 
search range within the model space, as well as the desired resolution (Stoffa y Sen, 1991). In 
this way, the range is defined for each parameter defining a couple of bounds a t and bi , that is 
to say, a t < $ < b t . The resolution is controlled with the discretization interval di , defined as 

= (b t - 

where Af,- is the number of possible values for the parameter during the process (Sambridge 
and Drijkoningen, 1992). The allowed models, e, defined by the set of parameters £ , are 
restricted to the domain of values 



s^a^jd, for y = 0,..., AT. (17) 
[0097] Initial population 

[0098] Also in block (1) of figure 7, from random combinations of the parameters, an 

initial population of models is constructed, whose dimension depends on the particular 
problem to be solved. Each combination translates into a set of integer indexes, defined by 
equation (17). Said integer values, which establish the particular value for each model 
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parameter, are later codified as binary chains called chromosomes. Said chains are made up 
of consecutive groups of bits, called genes, that represent the value of the different 
parameters £) . The models that form the population in later generations are created applying 
the three essential evolutionary mechanisms: selection, crossover and mutation. 
[0099] Forward problem and evaluation of the cost function 

[00100] The forward problem in block (2) of figure 7 consists in calculating the 

theoretical or synthetic response for each model during the iterative process. Said response is 
then compared in block (3) of figure 7 with the observations or data using some measure of 
similarity known as cost function. The criterion utilized in this work is the L2 norm already 
defined in equation (14), albeit there can be others. The norm chosen to evaluate the fit must 
take into account the form and complexity of the observed and calculated curves. 
[00101] Selection 

[00102] In block (4) of figure 7, starting from a population of Q individuals and their 

respective cost functions E(sk) (k = 1, ...,Q), a ('cumulative') probability of selection P(sk) is 
determined that will depend on its degree of fitness. 

[00103] One formula that can be used to determine the cumulative probability of 

selection is 

P(s k ) = P(s t _,) + *™~ E _ { l k) (18) 

z£ \ max avr s 

where E max , and E avr are maximum and average cost functions of the generation, respectively, 
and Q is the number of individuals in the population. Next, a biased roulette procedure 
(Goldberg, 1989) is used to select a new population of models. Q random numbers r k between 
cero and one are generated. If P(e k .\) <r k < P(Sk), then e* is selected to be part of the new 
population. This means that in the new population there may be some 'twin 1 models. 
Additionally, 'clones 1 of the best models can be added to the base population of Q models. 
Said 'clones' will not suffer crossover nor mutation, in order to ensure that these desirable 
models are not lost when going through those processes, which are of a random nature. 
[00104] Alternatively, the selection probability can be determined as (Sambridge y 

Drijkoningen, 1992) 

P(e k ) = a-bE{z k ) (19) 
which describes a linear probability distribution, and 
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P(s k ) = Ae- BE{Zk) (20) 

which corresponds to an exponential distribution. The values that are often given to the 
constants a, b, A y B are as follows 

b = Q' l (E - E Y l , a > bE 

is- \ max avr J ' — 



5 = {EX » ^ = 



mar 
-i-l 



where jE^™, 2s flvr and jE^are the maximum and average cost functions, and the standard 
deviation of all the misfits of the initial population, respectively. 
[00105] Stoffa and Sen (1991) propose a selection criterion based on an update 

probability. The criterion consists in comparing the misfit of each model from the current 
generation with that of a model from the previous generation, chosen at random. If the misfit 
of the current model is smaller, then it is saved. If not, a value P u is considered that 
establishes the probability of substituting the current model for the previous one. The 
procedure controls the influence of the misfit of previous generations on the current 
population. The value suggested by Stoffa and Sen (1991) for the substitution probability P u 
is 90%. 



[00106] Crossover 

[00107] The new models are produced in block (5) of figure 7 from a parent generation 

of Q models. Q/2 pairs are formed randomly. Each pair is potentially able to mate. In order to 
determine which pairs are to be crossed over, each one is assigned a random number between 
cero and one. If said number is smaller that the crossover probability P c , then the 
corresponding pair are crossed over. If not, both members of the pair are kept unchanged in 
the next generation. 

[00108] The crossover mechanism is based on choosing randomly the position of a 

gene for both chains. The chains are broken at that point to interchange information between 
them, as can be seen in figure 8, where is shown in (A) the partition of the chromosomes into 
four garnets, and in (B) the exchange and union of two gametes to form two cigots. The 
purpose of crossing over two different chains is to explore new regions of the solution 
domain, where the absolute minimum might lie. In the normal crossover process, the mated 
pairs have two children, and the population of models is automatically maintained in Q 
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individuals. Alternatively, it can be considered that the couples have only one child, in which 
case the population must be completed with randomly generated models until it reaches Q 
individuals. 



[00109] Mutation 

[00110] Mutation, like the sexual reproduction (or crossover), contributes to the 

genetic diversity of a population. Mutation allows the search to prosper when it is confined in 
the vicinity of a local minimum. The mutation is done in block (6) of figure 7 by inverting a 
bit chosen at random within the binary chain (chromosome). The number of models to which 
the process of mutation is applied, as in the crossover, depends on a parameter called 
probability of mutation, denoted by P m . This mechanism prevents the premature convergence 
of the method, when the population is too homogeneous and incapable to continue the 
evolutionary process. In figure 9 an arbitrary chain is represented and the mutation of a 
random bit is exemplified, the original chromosome is shown in (A) and the mutated one is 
shown in (B). 

[00111] An alternative to define the probability of mutation P m was proposed by 

Yamanaka and Ishida (1996). It consists in determining the level of homogeneity of the 
individuals in every generation through the calculation of an average variation coefficient^, 
for each parameter, using the formula 



P wU* 



(21) 



where p is the number of parameters, e g is the average of the z'-th parameter, and oi is the 
standard deviation. Below, P m is defined as a function of y 



P ini P^a y>0.\ 

0.1 para 0.02 <x< 0.1 (22) 
0.2 para y < 0.02 



where P ini is the initial probability of mutation. 

[00112] With mutation concludes the sequence of operations that define a genetic 

algorithm. Said sequence is repeated until some pre-established tolerance is satisfied. 
[001 13] As an example, one possible specific embodiment, albeit not the only one, of 

the method of genetic algorithms, is presented in figure 10. In block (1), the method starts 
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with a population of Q randomly generated models, each one representing a permittivity 
distribution. At the same time, also in block (1) the iteration counter i is initialized. Next, the 
model counter k is initialized and a cycle begins through which, for each one of the Q 
models, the forward problem is solved and the misfit function E(s k ) as well as the cumulative 
probability of selection P(em) are calculated in block (2). After that, selection of the models 
with lower E is performed in block (3) according to the value of P(Sk) and applying the 
process of biased roulette previously described. Then follows model crossover in block (4). 
In said block (4), Q/2 pairs of models are formed at random, which are then selectively 
crossed over based on the value of P c , producing only one child per pair. Also in block (4), 
the pairs not crossed over remain (both members) unchanged and the population is completed 
to Q individuals with new random models. Finally, in block (5) the average variation 
coefficient y and the mutation probability P m are calculated, and then model mutation is 
carried out according to the description already given in previous paragraphs. The whole 
previous procedure is repeated continuously, increasing i by one on every iteration. The 
algorithm ends when a pre-set number of iterations N is reached or if the misfit function is 
sufficiently small. 

[00114] Formulation of the Forward Problem 

[00115] The forward problem consists in calculating the mutual capacitances c,y , i j\ 

that result from the presence of a permittivity distribution e inside the sensor. Both methods, 
simulated annealing and genetic algorithms, require the repeated solution of the forward 
problem. Because of that, it is important to have a suitable method to solve said problem, that 
achieves a reasonable balance between accuracy (or precision) and speed. In the context of 
this invention, the forward problem can be solved using an optimized routine developed by 
the authors based on the finite-volume method, which will be described briefly. This routine 
is more efficient than those reported so far in the literature (Yang y Peng, 2003) for it is 
comparable in its precision with implementations based on the finite-element method using 
meshes with 9,000 triangular elements. The execution speed is higher than that of finite- 
element and finite-difference methods. The routine is written in Fortran 90 and is totally 
portable (it has been tested on PC-type computers, MS-Windows and Linux based PC 
clusters, SUN and ALPHA workstations, and CRAY supercomputers). The routine can be 
extended to the three-dimensional case without any major modifications and is easily 
parallelizable. 
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[00116] The forward problem is solved using the finite-volume method in a cylindrical 

configuration. In this way, the intermediate solutions in the center of the disc (which are a 
problem in the finite-difference method) are eliminated and the mesh refinement becomes 
more flexible compared with finite-element methods. Equation (7), repeated below, is solved 

V-s(x 9 y)V<f> k =0 

[00117] where e is the permittivity and <p k is the electrostatic potential generated when 

electrode k is the source (or excitation). The equation is subject to the boundary conditions 
(a) </> k = V volts on the source electrode and (b) (f> k = 0 on the detection electrodes and on the 
outer screen. 

[00118] Defining the radial and angular coordinates as r and 0, and using the finite- 

volume method, the discrete equation is formulated in conservative form for each cell ft^ as 



£ V-(fV^) JQ.. = 0 for i=\ 9 .. 9 N r and j=l 9 .. 9 N 0 



(23) 



where i and j refer to the discretization in r and 0 9 respectively, and N r and No are the number 
of sections into which the radius and the circumference are divided, respectively. 
[00119] Applying Gauss's theorem in polar coordinated, the discrete equations can be 

written as 



f sV<p k • dY iy = 0 



(24) 



where Yy is the boundary of the finite volume cell Q,y . The boundary r,y is defined by Y w and 
Y E along the radial coordinates, and by Y N and Y s along the angular coordinates. Equation 
(24) can be expressed as the sum of the fluxes through the faces Y N , F s , T E and Y w 



- e 



s d(f> k 
r 69 

d<f> k 



Ar 



80 



rA0 



(25) 



[00120] From (25), the term corresponding to the fluxes at cero radius vanishes and the 

problem is equivalent to solving the equations in the proximity of the center on triangles that 
have a vertex on the center. Then, the discrete system of equations for the forward problem is 
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well posed. The complete system is similar to a Laplacian system of equations, and a 
diagonal banded system that includes the periodic boundary conditions imposed by the 
problem geometry must be solved. The corresponding matrix is positive definite and non- 
symmetric, characteristics that can be exploited when selecting the specific optimum solution 
methods. 

[00121] Finally, the mutual capacitances are calculated by integrating the potential 

gradients along a curve surrounding the electrodes, according to equation (6), which is 
repeated below 




[00122] where e Q is the permittivity of free space (8.854 x 10" 12 farads per meter), T, is 

a closed curve surrounding electrode i, dl is a normal vector representing an element of the 
curve T/, dl is an element of length of that curve and (f> k is the electrostatic potential produced 
in the sensor when applying a voltage of V volts to electrode k (source or excitation) and 0 
volts to all others (detection electrodes). The integration is done using a trapezoidal rule and 
the potential gradients are calculated to the fourth order. 

[00123] During the procedure for reconstructing a permittivity image using simulated 

annealing, it is necessary to solve the forward problem and find the electric potential 
repeatedly for relatively similar successive permittivity distributions, while the method 
converges towards the final solution. Since the potential corresponding to said successive 
distributions changes relatively little, it is possible to accelerate the whole process if an 
iterative method is used to solve the forward problem, taking as the first guess for the 
potential the solution potential corresponding to the previous permittivity configuration. 
Because the initial guess for the potential will be quite close to the solution, said iterative 
method will converge in a few iterations, rapidly achieving an acceptable accuracy. 



[00124] EXAMPLARY EMBODIMENTS [[Ell 

[00125] With the purpose of evaluating the performance of the image reconstruction 

methods described previously, a set of synthetic ECT data were calculated using the forward 
problem routine. In order to emulate the main sources of uncertainty in the data produced by 
a measuring instrument operating in real working conditions (i.e., the random errors inherent 
to the process of measuring any physical quantity and the errors caused by the limited 
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precision of the measuring instrument), the calculation of the synthetic capacitances for the 
ideal model were evaluated with a numeric precision of the order of 10" 1 1 in the iterative 
method used in the forward problem for the calculation of the potential. To emulate the 
(systematic) imprecision associated with the ECT sensor, during the inversion process (i.e., 
the estimation of the electrical permittivity distribution inside the pipe) a considerably smaller 
precision (10~ 5 ) was used in the calculation of the potentials in the forward problem. 
[00126] The interpretation of potential-field data without any restriction can be of very 

little practical interest due to the great ambiguity present between the observations and the 
estimated solutions. The non-uniqueness in the potential-field problems arises mainly from 
two sources: the first is the inherent ambiguity caused by physics of the problem that permits 
the existence of many solutions that reproduce the anomaly in the potential field; the second 
results from the use of a finite number of data that are contaminated with errors and that may 
not contain enough information to construct a unique solution to the problem. The strategies 
that allow overcoming this non-uniqueness consist in the incorporation of sufficient a-priori 
information to constrain the resulting solutions to a region in the parameter space that is 
considered physically reasonable (Pilkington, 1997). In the particular case of electrical 
capacitance tomography, there is information about the typical values of the electrical 
permittivity for gas, oil and water. The knowledge of these values, as well as the significant 
contrast between the properties of gas and oil, allow the application of the inversion 
techniques discussed here to ECT data for two-phase flows contaminated with relative errors 
of up to 2% (this being the maximum error produced by the data acquisition systems 
normally employed in ECT). Additionally, if there is a precise statistical estimation of the 
data uncertainty, then it is possible to construct a model that considers the uncertainties in the 
data and in the parameter estimations, using the scheme proposed by several authors for the 
inversion of potential-field data (for example, Sen and Stoffa, 1995). 
[00127] The synthetic capacitances were calculated for three typical permittivity 

distributions using the forward problem subroutine based on finite volumes. A 12-electrode 
ECT sensor was simulated and the capacitance values for all possible electrode-pair 
combinations were calculated. Those were the simulated data. A two-component distribution 
was considered with a material with a low permittivity of 1 (gas) and another with a high 
permittivity of 2.5 (oil). Tests were carried out with noise-free data and with data 
contaminated with random errors of up to 1%. The algorithm is written in Fortran 90 and runs 
on a L7-gigahertz Pentium 4 computer with 512 megabytes of RAM memory. The tests were 
done using a mesh having 120><60 elements in order to reduce the inversion time (~ 30 
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minutes for 60,000 iterations), but the results are valid for larger dimensions. The validity of 
these simulations was corroborated with experimental laboratory results obtained using 
physical models and a real electrical capacitance tomography system. 
[00128] After a suitable parameterization, both simulated annealing and genetic 

algorithms yielded very similar results for the three cases studied. The quality of the 
reconstructed permittivity images depends mainly on the number of iterations of the method, 
as it happens in many other applications (for example, Ortiz Aleman et al., 1999, 2001, 2002, 
2003; Cruz-Atienza, 1999). In figure 1 1, the reconstructions for a simplified annular flow (A) 
are presented, after 30,000 (B) and 60,000 (C) calculations of the direct problem (i.e., 
iterations of the inversion method). In figures 12 and 13, similar illustrations are shown for a 
stratified flow and a bubbly flow, showing in (A) the ideal image, in (B) the reconstructed 
image after 30,000 iterations, and in (C) the reconstructed image after 60,000 iterations. In all 
cases shown in the figures, the simulated annealing method was used starting from an initial 
homogeneous permittivity distribution. These results show clearly that the reconstructed 
images closely resemble the 'true 1 ones. The accuracy of these reconstructions is considerably 
better than those reported so far in the literature (Yang and Peng, 2003). Albeit the methods 
of this invention do not require that the initial guess for the permittivity distribution be close 
to the solution in order to converge, it is possible to speed up the process a little if the image 
resulting from a simple direct method such as LPB is used as the first guess. 
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